ITEP-LAT/2010-07 



Schwinger-Dyson equations in large- N quantum field theories and nonlinear random 

processes 



O 

(N 

X> 
IX, 

(N 



P. V. Buividovich 1 ' 2 '* 

1 1TEP, 117218 Russia, Moscow, B. Cheremushkmskaya str. 25 
2 JINR, 14 1980 Russia, Moscow region, Dubna, Joliot- Curie str. 6 
(Dated: December 27, 2010) 

We propose a stochastic method for solving Schwinger-Dyson equations in large-iV quantum 
field theories. Expectation values of single-trace operators are sampled by stationary probability 
distributions of the so-called nonlinear random processes. The set of all histories of such processes 
corresponds to the set of all planar diagrams in the perturbative expansions of the expectation 
values of singlet operators. We illustrate the method on the examples of the matrix- valued scalar 
field theory and the Weingarten model of random planar surfaces on the lattice. For theories with 
compact field variables, such as sigma-models or non-Abelian lattice gauge theories, the method does 
not converge in the physically most interesting weak-coupling limit. In this case one can absorb the 
divergences into a self-consistent redefinition of expansion parameters. Stochastic solution of the 
self-consistency conditions can be implemented as a "memory" of the random process, so that some 
parameters of the process are estimated from its previous history. We illustrate this idea on the 
example of two-dimensional O (N) sigma-model. Extension to non-Abelian lattice gauge theories is 
discussed. 
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Introduction 



Modern lattice QCD simulations are mostly based on 
direct evaluation of the path integral of the theory. Such 
approach, while being very general and efficient for many 
applications, suffers from a number of problems, most no- 
table of which are the sign problem at finite chemical po- 
tential, critical slowing down at small quark masses and 
large finite- volume effects as well as small signal-to-noise 
ratio in the analysis of excited states. These problems are 
inherent to standard Monte-Carlo simulations and can- 
not be efficiently solved by simply increasing the compu- 
tation power, since the required computing time quickly 
increases (in the worst cases, exponentially) with the re- 
quired precision. Such situation makes it tempting to 
devise alternative simulation algorithms for non-Abelian 
lattice gauge theories. 

One of the efficient alternative numerical methods 
is the so-called Diagrammatic Monte-Carlo, a method 
based on stochastic summation of all the terms in the 
strong- or weak-coupling expansion of the observable of 
interest [1, 2]. Such a method in some cases allows one to 
reduce or avoid completely the sign problem in the origi- 
nal path integral, and does not suffer from finite- volume 
effects. Furthermore, one can construct algorithms which 
yield particular correlation functions in terms of prob- 
ability distributions of some random variables, which 
greatly facilitates the analysis of excited states [1, 2]. 
This is the idea of the "worm" algorithm by Prokof'cv 
and Svistunov [1], in which the probability distribution 
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of the positions x, y of "head" and the "tail" of the worm 
yields the two-point Green function G(x,y). Diagram- 
matic Monte-Carlo and the "worm" algorithm have been 
successfully applied to a number of statistical models 
with discrete symmetry groups such as the Ising model, 
the XY model and unitary Fermi gas and showed practi- 
cally no critical slowing down near quantum phase tran- 
sitions. 

However, application of such methods to lattice field 
theories with continuous field variables (such as two- 
dimensional O (N) and CP (N) sigma-models, Abelian 
gauge theories and the </> 4 theory) resulted so far in quite 
complicated and model-dependent algorithms [2] . A gen- 
eralization of such algorithms to SU (N) sigma-models or 
to non-Abelian gauge theories is still not found. These 
algorithms are in essence based on the strong-coupling 
expansion, and while their applicability is not limited 
by the strong-coupling regime, one can expect that algo- 
rithms based on the weak-coupling expansion might show 
better performance near the continuum limit. 

Typically, the weak-coupling expansion in such lattice 
theories is either quite complicated or non-convergent. 
Up to now, divergent behavior of the weak-coupling per- 
turbative expansions strongly limits the applicability of 
Diagrammatic Monte-Carlo to field theories with contin- 
uous field variables. In a recent paper [3] a method was 
proposed to construct convergent series which approxi- 
mate the non-analytic path integrals with desired preci- 
sion. This method, however, is difficult to generalize to 
physically interesting field theories such as non-Abelian 
lattice gauge theories. 

Another way to obtain convergent series while preserv- 
ing important physical properties of the theory is to sum 
over diagrams with certain topology only. This corre- 
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sponds to the largc-iV limit in quantum field theories 
and matrix models, that is, the limit of infinite dimen- 
sionality of an internal symmetry group, such as O (N) 
or SU(N). For such theories, each Feynman diagram 
acquires a factor N x , where x 1S the Euler character of 
this diagram [4]. In the limit N — > oo, the contribution 
of planar diagrams with x — 2 dominates, and the sum 
over all planar diagrams typically has a finite convergence 
radius [5, 6]. 

In this paper we describe a stochastic method for sum- 
ming over all planar diagrams in large- N quantum field 
theories. The method is based on stochastic solution 
of Schwinger-Dyson equations, so that the correlators of 
field variables are obtained as stationary probability dis- 
tributions of certain random variables. In this way we 
implement the idea of importance sampling, so that nu- 
merically small observables correspond to unlikely events. 
These probability distributions are sampled by the so- 
called nonlinear random processes. In contrast to con- 
ventional Markov chains, stationary probability distribu- 
tions of such random processes satisfy nonlinear equa- 
tions, and hence they can be called "nonlinear random 
processes" or "nonlinear Markov chains" in the termi- 
nology of [7, 8] . Factorization of single-trace operators in 
the large- N limit of quantum field theories corresponds 
to the phenomena of "chaos propagation" in random pro- 
cesses [9]. 

While in the diagrammatic Monte-Carlo and in the 
"worm" algorithm the diagrams are stored in computer 
memory as a whole and are updated in such a way that 
the detailed balance condition is satisfied at each step, 
the method described in this paper works only with ex- 
ternal lines. In contrast to the standard Metropolis algo- 
rithm, one should not know explicitly the weight of each 
diagram, and the transition probabilities do not satisfy 
any detailed balance condition. Unlike the quite popular 
"numerical functional methods" in continuum gauge the- 
ories (see [10] for a review), the proposed method does 
not require any truncation of the hierarchy of Schwinger- 
Dyson equations, and work only with singlet operators 
w.r.t. the internal symmetry group. Another distinct fea- 
ture is that the computational complexity of the method 
does not depend on N, while the standard Monte-Carlo, 
the functional methods and the "worm" algorithm all re- 
quire infinite computational resources in the limit N —> 
oo. This feature might be advantageous for numerical 
checks of the predictions of the holographic models which 
are dual to large- N quantum field theories [11]. 

In Section 1 we analyze the general structure of 
Schwinger-Dyson equations in large- N quantum field the- 
ories on the example of a scalar matrix- valued field the- 
ory. When large-iV factorization is taken into account, 
Schwinger-Dyson equations become nonlinear equations 
with infinitely many unknowns. In Section 2 we de- 
scribe nonlinear random processes of recursive type [7] 
which can be used to stochastically solve such equations. 
In Section 3 we apply such random processes to solve 
Schwinger-Dyson equations in several large-TV theories. 



In Subsection 3.1 we consider the scalar matrix- valued 
field theory, for which the perturbative expansion yields 
the conventional Feynman diagrams in momentum space. 
In Subsection 3.2 this solution is compared with the ex- 
act solution of the simplest quantum field theory in zero 
dimensions, that is, the Hcrmitian matrix model [6]. The 
convergence of such solution and the strength of the sign 
problem is discussed. In Subsection 3.3 we consider the 
Weingarten model [12, 13] and demonstrate how the pro- 
posed method can be used to simulate random surfaces 
on the hypercubic lattice. In this case, our method re- 
produces an ensemble of open, rather than closed, ran- 
dom surfaces, with critical behavior which is quite dif- 
ferent from those of the closed planar random surfaces. 
Since the structure of Schwinger-Dyson equations in the 
Weingarten model is similar to the loop equations in 
large- TV non-Abelian lattice gauge theories [14], study- 
ing this model might be helpful for further extensions of 
the present approach to non-Abelian gauge theories. 

While the method described in Section 2 works well 
for non-compact field variables, for field theories with 
compact field variables, such as nonlinear sigma-models 
or non-Abelian lattice gauge theories, a straightforward 
stochastic interpretation of Schwinger-Dyson equations 
is only possible in the strong coupling limit. In the 
weak-coupling limit one expects the field correlators to 
contain both the perturbative part in the coupling con- 
stant g as well as nonperturbative corrections of the form 
exp (— c/g 2 ) with some constant c. Moreover, perturba- 
tive expansion in powers of g typically results in asymp- 
totic series, and nonperturbative corrections appear as a 
result of resummation of such series [15]. 

In Section 4 we show how such nonperturbative correc- 
tions can be taken into account by a further relaxation 
of the Markov property of the random process. The ba- 
sic idea is to absorb the divergent part of the series into 
a self-consistent redefinition of the expansion parameter. 
These redefined parameters play the role of nonpertur- 
bative "condensates" [15, 16]. It turns out that the re- 
defined expansion parameters can be estimated with in- 
creasing precision from the previous history of the ran- 
dom process which solves the Schwinger-Dyson equa- 
tions, thus leading to the emergence of the "memory" 
of the random process. The approach of the redefined 
parameters to their self-consistent values is reminiscent 
somehow of the rcnormalization-group flow [10]. Such 
dependence on the previous history makes the random 
process essentially non-Markovian, so that the station- 
ary probability distribution also satisfies some nonlinear 
equation. 

We illustrate this idea on the example of O (N) sigma- 
model in two dimensions, which is equivalent to a bosonic 
random walk with a self-consistent mass. Random pro- 
cess which simulates this model has the "memory" but 
no "recursive" structure. Presumably, in order to sum 
up both perturbative and non-perturbative corrections 
which arise in non-Abelian lattice gauge theories or 
U (N) sigma-models, one should devise the "recursive" 
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nonlinear random process (which would sum up pertur- 
bative corrections) with memory (which would generate 
nonperturbativc quantities in a rcnormalization-group- 
like way). 

Finally, in the concluding Section we summarize the 
present work and discuss its extension to non-Abelian 
lattice gauge theories in the limit of large N. 

1. GENERAL STRUCTURE OF 
SCHWINGER-DYSON EQUATIONS FOR 
LARGE- iV QUANTUM FIELD THEORIES 



This theory is most convenient to illustrate the method 
described in this paper, since its perturbative expansion 
leads to conventional Fcynman diagrams in the momen- 
tum space. Since this theory should be somehow regu- 
larized, let us assume from the very beginning that the 
action (1) is defined on the Euclidean hypercubic D- 
dimensional lattice with total volume V in lattice units. 
Thus, the coordinates x take integer values and A is the 
lattice Laplacian (for definiteness, with periodic bound- 
ary conditions). 



In order to analyze the general structure of Schwinger- 
Dyson equations for large-iV quantum field theories, let 
us first consider the theory of a hermitian N x N matrix- 
valued field (ft (x) with the following Lagrangian: 

C [eft (a?)] = iVTr eft (x) (m 2 ~ A) eft (x) + ^Tr eft 4 (x) .(1) 



Schwinger-Dyson equations for a theory with the ac- 
tion (I) read [17]: 



(m 2 - Ai) G(x 1} x 2 ) = 5 (xi,x 2 ) + XG (x 1 ,xi,xi,x 2 ) , (2) 
(m 2 - Ai) G (xi, ...,x n ) = 6 (xi,x 2 ) G(x 3 , ...,x n ) + 



n-l 



+8 (x!,x n ) G (x 2 , ■ ■ .,a; n _i) + ^ S ( x ii x m) G(x 2 , . . .,x m -i) G(x m +i,- ..,x n ) + 
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h\G(xi,xi,xi,x 2 ,...,x n ), n>2, (3) 

where the single-trace correlators are G (x\, . . . , x n ) = ( ■h Tr (eft (x\) . . . eft (x n )) ), Ai is the Laplacian acting on x\ 
and we have already taken into account the factorization property in the limit N — > oo [4]: 

( ^Tr (eft ( Xl ) ... (x n )) ^Tr (eft ( Vl ) ... (y m )) ) = 

(4) 



( llY (eft ( Xl ) ...eft (x n )) ) ( lit (eft ( yi ) ...eft (y m )) ) + o(Jp 



These equations hold for any argument of the correlators, but the resulting system is redundant, and it is sufficient 
to consider only those Schwinger-Dyson equations which were obtained by the variation of the fields at x\. 

It is convenient now to go to the momentum representation, introducing the Green functions in momentum space 



G (A;i, ...,&„) = X} • • - 2 ex P ( * X ^™ ' x m ) G (xi, . . . , x n ). In order to keep all expressions as symmetric as possible, 
we do not separate the factor 6 k m I in G (k\, . . . , k n ) explicitly. This condition will be automatically satisfied 



by the nonlinear random process which we describe in Subsection 3.1. The equations (2), (3) in the moment 
representation are: 

G(k 1 ,k 2 ) = G (k 1 ) VS(k 1 + k 2 ) + 
+G (ki) X! 8 [ki - qi - qi - qs) G (qi,q 2 ,q 3 ,k 2 ) , (5) 
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G(ki, k n ) = G (fci) ^ S (h + k m ) V G(k 2 , fc m -i) G (k m+ i, k n ) + 

m— 3 

+G (k 1 )S(k 1 + k 2 ) VG(k 3 ,...,k n ) + G (h) 5(k 1 + k n ) VG(k 2 , fc„_i) + 



+G (k 1 ) — H k i-qi-qz-q3) G(qi,q 2 ,q 3 ,k 2 ,...,k„), (6) 
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where Go (k) = m 2 + 4 sin 2 (k^/2) is the free scalar propagator on the hypercubic lattice. All momenta are 
.A M . . / 

assumed to lie in the first Brillouin zone — tt < k^ < n and are added modulo 2tt. The structure of these equations 
is schematically illustrated on Fig. 1. where dashed blobs denote the Green functions G (ki, . . . , k n ) and empty blobs 
denote Go (k). 




FIG. 1: Schematic illustration of the structure of the 
Schwinger-Dyson equations (6). Dashed blobs denote the 
Green functions G (fci, . . . , k n ) and empty blobs denote the 
free propagator Go (k). 



Thus we have obtained an infinite system of 
quadratic functional equations for the set of functions 
G (fci, . . . , k n ) with n = 2,4,.... Such structure is com- 
mon for largc-iV quantum field theories: Schwinger- 
Dyson equations are quadratic equations for infinite set 
of unknown variables. In the case of scalar matrix field 
theory considered here, the unknown variables are the 
functions of the sequences of momenta {ki, . . . , k n } for 
any even n > 2. In the case of lattice gauge theories or 
string theories Schwinger-Dyson equations are most nat- 
urally formulated in terms of the Wilson loops, which are 
the functions defined on the discrete space of closed loops 
on the lattice [12-14]. In this case, the equations are also 
quadratic w.r.t. the Wilson loops. For O (N) sigma- 
model, Schwinger-Dyson equations are also quadratic 
equations which involve only the two-point function (see 
Section 4). 

Typically, systems of equations with infinitely many 
unknowns can be efficiently solved by stochastic meth- 
ods. It is advantageous to estimate the value of each un- 
known variable as a probability of observing some state 
of a random process. In this case the unknowns with nu- 
merically small values correspond to unlikely events, and 
the set of infinitely many unknown variables is automat- 
ically truncated to a set of unknowns with sufficiently 
large values. Such methods are well-known mainly in the 
context of kinetic equations [9] . Recently they were also 
discussed in the context of probabilistic programming [7] . 
In the next Section we describe a discrete-time, discrete- 
space method of such type which is in our opinion most 
suitable for solving the Schwinger-Dyson equations in the 
large- N limit. 



2. STOCHASTIC SOLUTION OF NONLINEAR 
EQUATIONS BY RANDOM PROCESSES OF 
RECURSIVE TYPE 

We consider nonlinear equations of the following form: 

w (x) = p c (x) + ^Pe (x\y) w (y) + 
y 

+ ^2 Pj (x\y l7 y 2 )w{y 1 )w{y 2 ) , (7) 

Vl,V2 

where x, y 7 yi, y 2 are the elements of some space X 
and ^2 with x € X denotes summation or integration 

X 

over all the elements of this space. We also assume that 
the functions p c (x), p e (x\y) and pj (x\yi, y 2 ) satisfy the 
inequalities 

\Pc (x) I + \ Pe (x\ yi ) I + \ Pj (x\ yil y 2 ) I < 1 (8) 

X 

for any yi, y 2 . 

We would like to find a stochastic process for which 
w (x) is proportional to the probability of the occurrence 
of the element x in some configuration space. Obviously, 
ordinary Markov process with configuration space X can- 
not solve such a problem, since stationary distributions 
of Markov processes obey linear equations. In order to 
solve the nonlinear equation (7), one can, for example, 
extend somehow the configuration space. Extensions of 
Markov processes with stationary probability distribu- 
tions which obey nonlinear equations have been consid- 
ered recently in [7, 8]. In this Section we concentrate 
on random processes similar to recursive Markov chains 
of [7]. The basic idea is that at any time one can leave 
the current chain and start a new one, then returning 
back to the old chain at some time. The initial state of a 
newly created chain depends on the states of older chains. 
Thus one has not a single Markov chain, but rather an 
infinite stack of chains. The random process which we 
describe below will be similar to these recursive Markov 
chains, but instead of referring to "recursion" we will ex- 
plicitly introduce the underlying stack structure. Here 
we first consider the equations (7) with the coefficients 
p c (x), p e (x\y) and pj {x\y\, y 2 ) being all positive, and in 
Appendix A we generalize to coefficients with arbitrary 
signs or complex phases. 

Consider an extended configuration space which con- 
sists of ordered sequences {x\, . . . , x n } for arbitrary n > 
1, with xi, . . . , x n € X . It is illustrative to interpret such 
configuration space as a stack of elements of the space X, 
so that x n is at the top of the stack. The desired random 



5 



process can be specified by the following prescriptions. 
At each discrete time step do one of the following: 

Create: With probability p c (x) create new element x G 
X and push it to the stack. 

Evolve: With probability p e (x\y) pop the element y 
from the stack and push the element x to the stack. 

Join: With probability pj (x\yi, y 2 ) consecutively pop 
two elements y\ , y 2 from the stack and push a single 
element x to the stack. 

Restart: With probability 1 

J2(Pc( x ) +Pe(x\yi) +p 3 (x\y 1 ,y 2 )), where y lt 

X 

2/2 are the two topmost elements in the stack, 
empty the stack and push a single element x £ X 
into it, with probability distribution proportional 
to p c (x). 

The last action is also the procedure used to initialize the 
random process. The "Evolve" action is just the evolu- 
tion of a single Markov chain at the top of the stack, with 
transition probabilities proportional to p e (x\y). The con- 
dition (8) and the positivity requirement ensures that 
p c (x), p e (x\y) and pj (x\yi, 1/2) can be interpreted as 
probabilities. 

Consider now an equation for the stationary probabil- 
ity distribution of such a Markov chain. It has a gen- 
eral form p(A) = ~^2,P(B — > A) p(B), where p(A) is 

B 

a stationary probability of the occurrence of a state A 
and P (B — > A) is the transition probability between the 
states B and A. Let W (xi, . . . ,x n ) be the stationary 
probability to find the elements x\,...,x n in the stack. 
This probability distribution function is obviously nor- 
malized to unity: 

00 

£X;...XV(zi,...s») = i. 

n— 1 x\ x n 

The equation for the stationary probability distribution 
in our case reads: 

W{x 1 )=N- 1 Pc (x x ) & + 5> e (a:i|l0 W (y) + 

v 

+ X] (zi|yi,2/2) W (#1,2/2) (9) 
W(xi, ...,x n )=p c (x n ) W (xi, . . . ,x„_i) + 

+ ^Pe {X n \y) W (Xl, . . .,X n -!,y) + 
V 

+ X p J ( x n\yi,V2) W{x 1 ,...,x n - 1 ,y 1 ,y 2 ),n > 1,(10) 

Vl,V2 



where 

-^Pe (x„\y) W(ari,...,x„_i,j/)- 
v 

- p i ( x n\yi>V2) W (x 1 ,...,x n -i,yi,y 2 ) 1(11) 

and jV c = ^Pc (x) . By a direct substitution one 

can check that there is a factorized solution for 
W (xi,. . .,£„): 

1¥ (ari, . . . , x„) = w (xi) w (x 2 ) ...w (x n ) , (12) 

where w (x) obeys exactly the equation (7) and wq (x) 
obeys the following inhomogencous linear equation: 

w (x) =Af~ 1 p c (x) ^ + ^p e (x\y) w (y) + 

y 

+ X] Pj( x \yi>V2) w (yi)w(y 2 ). (13) 

Thus, for any equation of the form (7) with positive 
coefficients which satisfy (8), there is a random process 
whose stationary distribution encodes the solution of this 
equation as in (12). The factorization of the stationary 
probability distribution of random processes with such 
an infinite configuration space is known as the "propa- 
gation of chaos" in random processes and was discovered 
for classical kinetic equations by McKean, Vlasov and 
Kac [9]. Comparing the equation (7) with Schwinger- 
Dyson equations (2), (3), (5) and (6), we conclude that 
this property corresponds to the factorization of singlc- 
trace operators in largc-iV quantum field theories. It is 
interesting that time reversal of the random process de- 
scribed above leads to the so-called branching random 
process [9], which has quite different properties. This is 
due to the fact that for such random processes there is 
no detailed balance condition, and hence no time reversal 
symmetry. We do not consider here a subtle mathemati- 
cal question of the existence of solutions to equation (7), 
since in our case it is ensured by the physical applications 
of this equation. 

Finally, let us describe a practical procedure for finding 
w (x) by simulating the random process described above. 
By standard statistical methods, one should sample the 
probability distribution p (x n ) of the topmost element in 
the stack (provided there is more than one element in it, 
otherwise we estimate wo (x) rather than w (x), see (12)). 
From (12), we get p (x n ) = s~ 1 w (x n ), with s = ^w (x). 

X 

It should be stressed that w (x) is not normalized to unity, 
but rather satisfies the inequality ^w (x) = s < 1. The 

X 

value of the normalization constant s can be also easily 
found numerically, since the probability to find n ele- 
ments in the stack decreases as s" for n > 1. 
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3. STOCHASTIC SOLUTION OF 
SCHWINGER-DYSON EQUATIONS BY 
RECURSIVE RANDOM PROCESSES 

3.1. Scalar matrix field theory 

After presenting the general method in Section 2, we 
are ready to describe a stochastic numerical solution of 
Schwingcr-Dyson equations (5), (6). For simplicity, let us 
assume that the coupling constant A in (1) is negative. 
This allows us to apply directly the results of Section 2, 
where all the coefficients in (7) are assumed to be posi- 
tive. In the case of positive A, additional sign variables 
for each sequence of momenta can be easily introduced 
following Appendix A. This will be done in the next Sub- 
section for the Hermitian matrix model. Note that while 
at finite N the theory with negative coupling constant is 
not defined and the correlators are non-analytic in A [3], 
in the leading order in N perturbative series converge 



Comparing the Schwingcr-Dyson equations (15), (16) 
with the general equation (7), we arrive at the random 
process which stochastically solves these equations. This 
random process is specified by the following probabilistic 
choice of actions at each discrete time step: 

Create: With probability Go (fc) (A/V) -1 push a new se- 
quence of momenta {fc, — fc} to the stack. 

Add: With probability Go (k) c~ 2 /V modify the top- 
most sequence of momenta {fci, . . . , k n } in the stack 
by adding a pair of momenta {k, —k} either as 
{k, fa, . . . , k n , fc} or {fc, fc, fa, . . . , kn}. 

Create vertex: With probability 
|A| Go (qi + q-2 + 53) c 2 replace the topmost se- 
quence {qi, q 2 , <?3, k 2 , • • • , k n } in the stack by 
{qi + q 2 + ?3j fa, • • • j kn}- This action can only be 
performed if the topmost sequence contains more 
than two elements. 



even when the coupling is negative, but not exceeding 
some critical value [6]. Correspondingly, in the planar 
approximation the correlators are analytic in A. 

The space X in (7) should be the space of ordered 
sequences (of any size) of momenta {fci, . . . , k n } , corre- 
spondingly, the extended configuration space is a stack 
which contains such sequences. It is convenient also to in- 
troduce two normalization constants M and c, so that the 
functions w (fa, . . . k n ) which will be estimated stochas- 
tically arc defined as 

G(fa,...,k n ) = J\fV n a™" 2 w (fa, . . . ,k n ) , (14) 

where V is again the total volume of space. The constant 
c can be thought of as the renormalization constant for 
the one-particle wave functions, and M - as the overall 
wavefunction normalization. 

In terms of the functions w (fci, . . . , fc„) the Schwingcr- 
Dyson equations (5) and (6) read: 



(15) 



(16) 



Join: With probability Go (fc) N c~ A /V pop the two se- 
quences {fci, . . . , fc„}, {q±, . . . , q m } from the stack 
(provided there are more than two elements in 
it) and join them into a single sequence as 
{fa, . . . ,k n ,k,qi, . . . ,q n ,—k}. Push the result to 
the stack. 

Restart: Otherwise restart with a stack containing a se- 
quence {fc, — fc}, fc being distributed with the prob- 
ability proportional to Go (fc) 

Since the momenta are always added to the stack in pairs 
which sum up to zero, for all sequences in the stack the 
total sum of all momenta in the sequence is always zero. 
The V~ x factors in (15), (16) ensure that the probabil- 
ity distributions of the newly created momenta can be 
normalized to unity. 

Let us check whether the inequalities (8) are satis- 
fied for such a process, that is, whether the total prob- 



w 



(fc 1 ,fc 2 ) = Go(fci)A ^ (fcl + fc2) + 



V 



G (fci) Ac 2 ^2 <K fc i - ?i - ?2 - 93) w(qi,q2,q3,k 2 ) , 

w(k 3 ,...,k n ) + 
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w(fa,...,k n ) = G (fci) c 



„ 2 S (fa + fc 2 ) 

V 



+G (fci)c — w (k 2 , ■ ■ ■ , k n -i) + 

+ Gq (fa)Nc~ 4 ^ ~~^T7 ~ W ( fc 2, ■ • ■ , fcm-l) W (fc m +l, ■ • • , fc„) 



m— 3 



V 



-G (fci) Ac 2 ^ Hfa -1i -92-93) w (qi, q 2l q3,k 2 , ...,k n ). 

I 
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ability of all possible actions docs not exceed unity. 
For the free propagator Go (k) one has the inequali- 
ties G (k) < l/m 2 and £ G (k) < V/m 2 . The total 

k 

probability of all possible actions can be then estimated 
as (Af- 1 + c- 2 + \\\c 2 +Afc- 4 ) /m 2 . Clearly, for suffi- 
ciently small |A| this estimate can be always made smaller 
than unity by increasing c and Af. In Subsections 3.2 and 
3.3 we will analyze such bounds on coupling constants 
in more details for Hcrmitian matrix model and for the 
Weingarten model. 

Since the constructed process involves no permuta- 
tions, one can trace the history of each momenta in the 
stack - from creation to joining into a vertex or a "Restart 
operation" . By drawing all the momenta in stack as 
points on the vertical lines of some two-dimensional grid 
and connecting the corresponding points along the hor- 
izontal lines, all planar diagrams of the theory (with an 
arbitrary number of external lines) can be obtained. Note 
also that the number of vertices in planar diagrams drawn 
by this random process cannot exceed the number of time 
steps from the previous "Restart" action. Thus in or- 
der to maximize the mean order of diagrams which are 
summed up in some fixed number of time steps, it is 
advantageous to maximally reduce the rate of "Restart" 
events, that is, to saturate the inequalities (8). 

One could also try to devise a random process which 
would solve the Schwinger-Dyson equations (2), (3) di- 
rectly in physical space-time, rather than in the momen- 
tum space. The configuration space of such a process 
would be the stack of sequences of points {xi, . . . ,i„}. 
As compared to the algorithm in the momentum space, 
there would be an additional choice of moving the last 
point x n in the topmost sequence to adjacent lattice sites, 
with the probability proportional to the hopping parame- 
ter k = (2D + m 2 ) . This would correspond to drawing 
the worldlines of virtual and real particles by bosonic ran- 
dom walks. Interestingly, such worldlines can be mapped 
onto the string worldshccts in simplicial string theory 
[18]. As well, the creation of a new interaction vertex 
would only be possible if three such random walks would 
intersect in one point. However, this is an unlikely event, 
with probability going to zero in the continuum limit. 
Thus, solving the Schwinger-Dyson equations directly in 
the coordinate representation would lead to a less effi- 
cient numerical algorithm. 

Note that for the theory (1) at finite N the 
Schwinger-Dyson equations are linear equations, which 
are, however, defined on much larger functional 
space: the set of unknown functions includes also 
expectation values of multi-trace operators, such as 
(Tr (<j>( Xl )...<t>(x n )) Tr (0 ( Vl ) . . . (y m ))) . One can 
try to solve these linear equations by interpreting them 
as the equations for the stationary probability distri- 
bution of a Markov process. The configuration space 
of such a process should be a space of sequences of 
the form {{xi , . . . , x n } , . . . , . . . , y m }}, thus encod- 
ing the expectation values of all multi-trace operators. 



However, such a straightforward procedure leads to non- 
normalizable transition probabilities, indicating that the 
series which one tries to sum up are divergent. Only 
when the terms subleading in 1 /N are omitted from the 
Schwinger-Dyson equations, they can be interpreted as 
stochastic equations. At the same time, we obtain the 
Markov process on the extended configuration space de- 
scribed in Section 2, which we interpret as the stack of 
sequences. The property of the "propagation of chaos" 
[9] ensures large- N factorization of single-trace operators 
(see equation (12)). We are thus led to the random pro- 
cess of recursive type [7] . 



3.2. Hermitian matrix model 

To check the considerations of the previous Subsection, 
let us consider the theory (1) in zero dimensions, that is, 
the hermitian matrix model with the following partition 
function: 

Z (A) = j \ J -/o, exp ^-A72Tr</> 2 + ^ Tr^ .(17) 

The Green functions now depend only on one integer 
n: G n = G{n) = (^Tr0 2 ™). The Schwinger-Dyson 
equations (2), (3) also take a very simple form: 

G x = 1 + AG 2 , 

n-2 

G n = 2G n — i + ^ ' G m G n — m — i + AG n+ i, n> 1.(18) 

m— 1 

Here we will assume that the coupling constant A can 
be both positive and negative, in order to illustrate the 
method described in Appendix A. Let us again define the 
"renormalized" Green functions w n as G n = Afc"^ 1 w n . 
In the case of arbitrary sign of A, the configuration space 
of the random process should be the stack which con- 
tains integer positive numbers and additional sign vari- 
ables. Following Appendix A, we introduce the variables 
and w n ' which are proportional to the probabili- 
ties to find the elements {n, +} or {n, — } at the top of 
the stack (provided the stack contains more than one el- 
ement). Then w n = — w n . We thus arrive at 
the following random process for stochastic evaluation of 
w n ■ At each discrete time step one performs at random 
one of the following actions: 

• With probability Af^ 1 add new element {1, +} to 
the stack. 

• With probability 2c _1 increase the topmost ele- 
ment in the stack by 1 and do not change its sign. 

• With probability | A| c decrease the topmost element 
in the stack by 1 (if it is greater than one) and 
multiply its sign by the sign of A. 



• With probability N c~ 2 pop the two elements 
{n, Si} and {m, s 2 } from the stack (provided there 
are more than two elements) and push the element 
{n + in + 1, S1S2} to the stack. 

• Otherwise empty the stack and push into it a single 
element {1, +}. 

Note that for positive A elements with the minus sign are 
not generated, so that Wn = and the random process 
automatically reduces to the one described in Section 2. 
The inequalities (8) read now: 

A/"" 1 +2C- 1 + \X\c + Mc- 2 < 1, Af>0, c>0(19) 

As discussed in Subsection 3.1, in order to increase the 
efficiency of the algorithm it is advantageous to saturate 
this inequality. It is easy to see that at the same time 
we saturate the upper bound on the absolute value of 
the coupling constant A. Maximizing this upper bound 
with respect to Af and c, we see that |A| cannot exceed 
the value A = 1/16 = 3/4 A c , where A c = 1/12 is the 
convergence radius of the planar perturbative expansion 
which can be found from the exact solution of the matrix 
model (17) [6]. Thus the described random process cov- 
ers only some finite subrange of coupling constants for 
which the model (17) is defined. It is easy to understand 
the origin of this limitation: in fact the random process 
described above simulates an ensemble of diagrams with 
an arbitrary number of external legs, with the weight of 
each diagram being proportional to A^" , where N v is the 
number of vertices. The number of open diagrams with 
a given number of vertices is obviously larger than the 
number of closed diagrams, hence the sums over open 
diagrams have smaller convergence radius. 



and express Af as a function of A: 



c 

(3 




FIG. 2: Green functions G„ (A) in (18) versus the coupling 
constant A for n = 1, . . . , 5, obtained after N = 10 6 discrete 
time steps of the described algorithm at fixed c = 8 and with 
Af given by "Branch 1" of (20). The error bars are smaller 
than the symbols on the plot. Exact results of [6] are plotted 
with solid lines. 



Af 



c 2 - 2c - |A|c 3 ± ^c 3 (c|A|-l)(4-c + c 2 |A|) 



(20) 



We call the solution with the minus sign in front of 
the square root "Branch 1" and the other solution 
"Branch 2" . 
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FIG. 3: Autocorrelation time of the random process described 
in Subsection 3.2 as the function of the coupling constant A 
at fixed c = 8 and for different choices of Af in (20). 

On Fig. 2 we plot the Green functions G n (A) evaluated 
using the described random process as functions of A up 
to n = 5. These results were obtained after N = 10 6 
discrete time steps at fixed c = 8 and with Af given by 
the "Branch 1" of (20). The error bars are smaller than 
the symbols on the plot. Solid lines are the exact results 
for G n (A) in the planar approximation, obtained using 
the saddle point method [6]. 
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FIG. 4: Mean stack size of the random process described in 
Subsection 3.2 as the function of the coupling constant A at 
fixed c = 8 and for different choices of Af in (20). 



For |A| < A, there is a continuous set of Af, c which 
saturate the inequality (19). One can, for example, fix c 



Autocorrelation time and mean stack size for the de- 
scribed random process are plotted on Figs. 3 and 4, 
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respectively, as the functions of the coupling constant 
A. The observable used to define the autocorrelation 
time was the sum of all numbers in the stack. First, 
we note that "Branch 1" is more advantageous for sim- 
ulations, since with larger mean stack size one can gain 
more statistics. However, in this case the autocorrelation 
time is also larger. Interestingly, for this branch both the 
autocorrelation time and the mean stack size have max- 
imum near A = rather than near the "critical point" of 
the random process A. For "Branch 2", these quantities 
increase slowly towards A. 



1 

0.8 



0.6 
0.4 
0.2 





I Br. 1, n = 1 
Br. 1, n = 3 
Br. 1 , n = 5 
Br. 2, n = 1 
Br. 2, n = 3 
Br. 2, n = 5 



-0.8 



-0.6 



-0.4 



-0.2 



X/\ r 



FIG. 5: The quantity in (21) for n = 1,3, 5, as the function 
of the coupling constant A at fixed c = 8 and for different 
choices of J\f in (20). "Br. 1,2" is for "Branch 1,2". 

In order to characterize the strength of the sign prob- 
lem, we consider the quantity 



,(+) 



.,(-) 



,(+) 



..(-) 



(21) 



= 1 if the random process generates only elements 
with the plus sign and = if the numbers of pluses and 
minuses exactly cancel. In practice, it is advantageous 
to have as large as possible, so that the difference 

u4 — Wn ' can be estimated with maximal precision. 

are plotted on Fig. 5 as the function of A for n = 
1, 3, 5. For A < 0, £ s (A) decreases with A and n. The sign 
cancelation is thus moderate for n = 1 (£f (—A) ps 0.75) 
and becomes more and more important for higher-order 
correlators - £| at A = —A is close to zero. It is interesting 
that are almost equal for two different choices of Af in 
(20). 



Thus there are no indications of severe critical slowing 
down in the whole range of possible coupling constants 
— A < A < A. The sign problem is also moderate for low- 
order correlators, but becomes more severe for higher- 
order correlators. It could be extremely interesting to 
extend the applicability of the described random process 
up to A = A c while preserving these attractive features 
of the algorithm. 




P 

b 



FIG. 6: Schematic illustration of the structure of the loop 
equations (24) in the Weingarten model (22). These equa- 
tions should hold for any link (marked by a thick line) which 
belongs to the loop. 



3.3. Random planar surfaces: the Weingarten 
model 



Weingarten model [12, 13] is a lattice field theory which 
in the large- N limit reproduces the sum over all closed 
surfaces with genus one on the hypercubic lattice. The 
action for each surface is proportional to its area, thus 
the model can be considered as a lattice regularization of 
bosonic strings with Nambu-Goto action. Although this 
model docs not have a nontrivial continuum limit for any 
space dimensionality [19], the structure of the functional 
integral and of the Schwinger-Dyson equations in this 
model are similar to those in large- N non-Abelian lattice 
gauge theory, and the analysis of this model might be 
helpful for the extension of the approach described here 
to non-Abelian gauge theories. In order to derive the 
Schwinger-Dyson equations, it is convenient to consider 
the reduced Weingarten model [13], which in the large-N 
limit is equivalent to the original model, similarly to the 
Eguchi-Kawai model for non-Abelian lattice gauge the- 
ory. It can be shown that in contrast to reduced lattice 
gauge theories, for the reduced Weingarten model addi- 
tional twisting is not necessary [20]. 

The reduced model is defined by an integral over com- 
plex N x N matrices J7„ = U_ u with fi = 1, . . . , D: 




If one treats the second term in the exponent in (22) as a perturbation and expands Z (/?) in powers of j3, the 
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resulting sum over planar diagrams is equivalent to the 
sum over all possible closed surfaces of genus one on the 
lattice with the weight /3' s ', where \S\ is the area of each 
surface. 

A basic observable in this model is the sum over all 
planar surfaces which are bounded by some closed loop 
C. The loop C can be uniquely specified by a sequence 
{fii, . . . , (J-n}, where /i's take the values ±1, . . . , ±D. In 
order to reconstruct the loop C from the sequence, one 
should start from an arbitrary point on a hypercubical 
lattice and move along one link in the direction fii, for- 
ward if /xi is positive and backward if fii is negative. 
From this new position one should similarly move in the 
direction fi 2 , and so on. From the diagrammatic expan- 



sion one can see that such a sum over surfaces is given 
by the following correlator: 

W(C) = W(m,... > n n ) = (^Tr (U^... C7 Mn ) ), (23) 

where one takes the conjugate variable [/J , if (j,a is neg- 
ative. This observable is similar to Wilson loop in lattice 
gauge theory but, unlike the Wilson loop, it does not 
have a "zigzag symmetry" [11]: passing a link forward 
and immediately backward changes the value of the Wil- 
son loop. In the large-TV limit the single-loop observ- 
ables factorize, which allows us to obtain a closed set of 
Schwinger-Dyson equations for W ifa, . . . , fi n ) [12, 13]: 



W(iii,..., /in) = 6(jn, -fi 2 ) W ifiz, ■ • • , /in) + S ini,-/j, n ) W i(l2, ■ ■ ■ ,fin-l) + 

n-1 

+ / ] W jll2,..., flA-l) W illA+1, ■ • ■ , fin) 6 (Ml! -fiA) + 
A=3 

+0 ^2 W ifi,fii,-fi,fi2,-.., fin), n>2. (24) 

IH^lMil 

These equations should hold for any lattice link fik belonging to the loop C, but the resulting system of equations is 
redundant, and it is sufficient to consider only one link [i\ on the loop. The equations (24) are schematically illustrated 
on Fig. 6, where the link /zi is marked by a thick line. 

We see that the equations (24) again take the form similar to (7). Let us now define the "renormalized" observable 
w ifii, . . . , fi n ) by rescaling W (/ii, . . . , /x„) by the factors M and q as W ifii, . . . , fi n ) = Nq n w ifii, . . . , fx n ). One can 
interpret the factor q" as the mass attached to the boundaries of random surfaces, somewhat like the bare quark mass 
in QCD. The equations (24) then take the following form: 

wifii,fi2) = (Nq 2 ) 1 8im,-ii 2 ) + Pq 2 ^ w ifi, fii,-fi, fi 2 ) 

W ifli, . . . , fl n ) = (f 2 5 ifll, ~fi 2 ) W ifls, ...,fl n ) + q~ 2 6 ifli, -fin) W (/i 2 , ■ • ■ , Mn-l) + 

n-1 

+Nq~ 2 ^ w (^2, ■ • ■ , flA-l) w if-A+1, ...,fln)6 (ill, ~ll A ) + 

+Pq 2 ^2 u, (/ i '^i5 - A t 5/ i 25---)Mn)) n>2. (25) 
ImI^ImiI 

I 



Let us now devise a random process of the type de- 
scribed in Section 2, which solves stochastically these 
equations. The configuration space is now a stack which 
contains closed loops, that is, sequences of indices h = 
±1, . . . , ±D. The desired random process is defined by 
the following possible actions at each discrete time step: 

Create a new loop: With probability 2DAf~ 1 q~ 2 cre- 
ate a new elementary loop C = {/j, —ft}, where 
H = ±1, . . . , ±D is random (cither positive or neg- 
ative) . 



Join loops: With probability 2DAfq 2 pop the two 
loops Ci = {fii, . . . ,fi n }, C 2 = {vi,...,v m } 
from the stack and form a new loop C by join- 
ing the loops Ci, C 2 with a link in the random 
direction /i (cither positive or negative): C = 
{(Mi, . . . , fi n , fx, vi, . . . , i/ m , — fj,}. This action can 
only be performed if there arc more than two loops 
in the stack. 

Flatten loop: If the three links in the end of the se- 
quence on the top of the stack form a boundary of 
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FIG. 7: The observables Wi x0 = W(fi,-fi) (1x0 loop) 
and Wixi = W (fi, u, —fx, —v) (lxl loop) in the Weingarten 
model as a functions of the coupling constant p for different 
dimensions D. Solid line corresponds to the first two terms in 
the perturbative expansion: W (fx, — /x) = 1 + 2 (D — 1) f3 2 + 
O (/3 4 ) , W {/J,, v, -fx, -v) = /3 + 8 (D - 1) p 3 + O (/3 5 ) . The 
data were obtained after 10 7 iterations of the algorithm de- 
scribed above. 



the plaquette, that is, if the topmost loop has the 
form C = {fii, . . . , fx n , fx,v, —fx} for some (X and v, 
replace these three links by a single link in the direc- 
tion v with probability fiq 2 : C = {fxi, . . . , fx n , u}. 

Append to loop: With probability ADq~ 2 append a 
pair {fi, — fi}, where fx is random (either pos- 
itive or negative), to the topmost sequence 
{fx\, . . . , fin} in the stack as {fx, —fx, fix, . . . , fx n } or 
{fx, fix, ... , fi n , —fx}. The probabilities of these two 
choices are equal. 

Restart: Otherwise start with a stack containing an el- 
ementary random loop C = {fx, —fx}, where ix = 
±1, . . . , ±D is chosen randomly. 

Again assuming that the sum of the probabilities of 
all possible actions is equal to one and the probability 
of "Restart" events is minimized, we obtain an equation 
relating /3, Af and q: 



Pq 2 + 2Dq- 2 (AT + AT' 1 + 2) = 1 



(26) 



Maximization with respect to Af yields the relation be- 
tween q and /?: 



<1 



1 ± VI - 32D/3 



(27) 



We call the solution with the minus sign in front of the 
square root "Branch 1" and the other solution "Branch 
2". 

The value of f3 in (27) cannot exceed the critical value 
(3 (D) = 1/ (321?) . As we have already seen on the exam- 
ple of the matrix model, this critical value does not neces- 
sarily coincide with the true critical point /3 C (D) at which 



FIG. 8: Two simple configurations of branched polymers on 
the lattice. The configurations on the left and on the right 
count as different configurations. 



the sum over planar surfaces diverges. Indeed, f3 (D) does 
not exceed the lower bound /3 C (D) > (24 (D — 1)) ob- 
tained in [13], and is significantly lower than the critical 
values obtained numerically in [20]. In fact, for P = pn 
all the observables are still dominated by the lowest-order 
perturbative contributions. Expectation values of the ob- 
servables Wixo = W (fx, —fx) (1x0 loop) and Wi x i = 
W {fx,u,—fx,—u) (lxl loop), which were obtained af- 
ter 10 7 iterations of the random process described above 
(with q given by "Branch 1" of (27)), are plotted on Fig. 
7 as the functions of the coupling constant p. Solid 
line corresponds to the first two terms in the pertur- 
bative expansion: VF ix o = 1 + 2 (D - 1) /3 2 + O (/3 4 ), 
Wixi = P + 8(D-1)P 3 + 0(P 5 ). Within statistical er- 
rors, one sees only the lowest-order perturbative contri- 
butions. It should be stressed that the proposed random 
process implements stochastic summation of diagrams of 
all orders, but due to the smallness of /3 < (3 (D), a very 
large computational time is required to see the contribu- 
tions of higher-order terms. 

Note that when the coupling constant p tends to zero 
(that is, the "bare string tension" of the random surfaces 
tends to infinity) and q lies between the two solutions of 
(27), the above random process describes just the growth 
of "branched polymers" , whose branches are bosonic ran- 
dom walks and hence correspond to particles rather than 
"strings" . These branches consist of loops in which ev- 
ery lattice link is passed twice and which hence sweep 
out zero area. 

Taking the limit P — > in (27), we find that the min- 
imal value of q is q = %/8-D. In order to understand 
this critical value we first note that in the limit P — >• 
the observables W {fx\, . . . ,fx n ) are all equal to one if the 
links fx\, . . . , fx n form a loop which sweeps out zero area 
and zero otherwise. The probabilities to encounter such 
loops in the described random process is hence propor- 
tional to w (fix, . . . , fx n ) ~ q~ n . Simple examples of such 
loops, which can be also thought of as the random tree- 
like graphs on the lattice, are shown on Fig. 8. Now 
imagine adding to some loop k links stemming from some 
lattice site. Since the loop includes each link twice, the 
probability decreases by q~ 2k . The number of possible 
configurations of k links is (2 x 2D) k since each of k 
links can point along any of D directions both forward 
and backward. An additional factor of 2 appears since 
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the zero-area loops pass twice through each point and 
the new links can be inserted between the links point- 
ing either forward or backward (for example, compare 
the configurations on the left and on the right of Fig. 
8). Finally, one can add any number k = 1,2,. ..,00 of 
branches to any point belonging to the branched poly- 
mer. At the criticality, adding any number of random 
links to some configuration should not change its over- 
all weight. Therefore, the change of the weight due to 
the added links times the number of ways to add them 
should be equal to unity. We are thus led to the following 
equation for q: 



+00 



ADq- 



k=l 



1 - ADq- 



(28) 



or 8Dq~ 2 = 1. Thus, in the limit j3 — > we indeed re- 
produce branched polymers with the correct critical be- 
havior. 

At nonzero j3, deviation from trivial branched poly- 
mer configurations can be characterized by the rate of 
the "Flatten loop" events. Indeed, since the probability 
of n such events is proportional to (3 n , such sequence of 
events corresponds to a random surface (which is in gen- 
eral open) consisting of n lattice plaqucttcs, plus some 
number of random trees. One can therefore think of the 
described random process as of the process of drawing 
random loops which sweep out random planar surfaces. 
The average rate of "flattening" events is plotted on Fig. 
9 on the right as a function of the coupling constant (3 
for different dimensions D and for different choices of q 
in (27). One can see that in the whole range of coupling 
constants < (3 < (3 (D) the rate of flattening events 
is numerically very small. On the other hand, the num- 
ber of links in the loops, as well as the number of loops 
stored in the stack, are quite large. Mean stack size and 
mean length of the topmost loop in the stack are plot- 
ted on Fig. 10 as a function of the coupling constant /3 
for different dimensions D and for different choices of q 
in (27). One can conclude therefore that the "branched 
polymers" actually dominate in the properties of the de- 
scribed random process. Critical behavior of these ran- 
dom trees is universal for any dimension D, that is why 
such observables as the mean stack size or the mean loop 
length, which are mainly sensitive to the length of loops 
rather than to the area of random surfaces, practically 
do not depend on space dimensionality. 

While the closed planar surfaces in the vicinity of the 
true critical point /3 C (D) of the Weingarten model are 
also dominated by "branched polymers" [19], in our en- 
semble of open random surfaces this dominance can be 
thought of as the manifestation of the tachyonic insta- 
bility of open, rather than closed, strings. The fact that 
the critical coupling /3 (D) in our case is smaller than the 
true critical point (3 C (D) can be explained by the fact 
that the number of open surfaces with a given area is 
obviously larger than the number of closed surfaces with 
the same area. The true critical coupling constant /3 C (D) 
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FIG. 9: Mean rate of "flattening" events for the random pro- 
cess solving the loop equations in the Weingarten model as a 
function of the coupling constant /3 at different dimensions D 
and for different choices of q in (27). "Br. 1,2" is for "Branch 
1, 2". 



can be quite easily found by a very simple re-weighting 
procedure, which will be described in details in a separate 
publication. 



4. RESUMMATION OF DIVERGENT SERIES 
AND RANDOM PROCESSES WITH MEMORY 

In Section 3 we have described a stochastic method 
for the solution of the Schwinger-Dyson equations for 
theories with noncompact variables. This method works 
only at small coupling constants and implements stochas- 
tic summation of perturbative series. For theories with 
compact variables, such as nonlinear er-models or lat- 
tice gauge theories, the structure of the Schwinger-Dyson 
equations is such that the described method can be 
straightforwardly applied only in the strong coupling 
regime, where one can stochastically sum all terms in 
the strong-coupling expansion. However, the continuum 
limit of such theories typically corresponds to the weak- 
coupling limit. An additional complication is that for 
physically interesting theories the observables cannot be 
expressed as convergent power series in the small cou- 
pling constant g, but rather contain non-analytic part 
which is typically of the form exp (— c/g 2 ) with some con- 
stant c [15]. 

In this Section we point out one possible way to deal 
with this problem. The basic idea is to absorb non- 
perturbative corrections into some self-consistent redefi- 
nition of the expansion parameter [15, 16]. Recently, a 
similar resummation method was also considered in [21]. 
Solving the self-consistency condition leads to the con- 
cept of a nonlinear random process with memory [8], in 
which all previous history of the process is used to esti- 
mate the value of the self-consistent expansion parame- 
ter. 
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FIG. 10: Mean stack size (on the left) and mean length of the topmost loop in the stack (on the right) for the random process 
solving the loop equations in the Weingarten model as a function of the coupling constant /3 at different dimensions D and for 
different choices of q in (27). 



Let us illustrate this idea on the simplest example of 
O (N) sigma model in the limit of large N. The model 
is defined by the following path integral over unit TV- 
component vectors n (x) living on the sites of the D- 
dimensional hypercubic lattice: 



Z = J T>n{x) cxp 

\n(x)\=l 



N 



n {x) -n(y) 



A ^ 

<xy> 



(29) 



where summation goes over all neighboring lattice sites. 
Despite its simplicity, this model in D = 2 dimen- 
sions is asymptotically free and has a mass gap which 
depends nonpcrturbatively on the coupling constant 
A. Schwingcr-Dyson equations in this theory can be 
written in terms of the two-point function £(x,y) = 
(n(x) ■ n(y)), £(x) =£(x,0) as: 

e(z) = ^£«( x± ^) (±e„)) + S (x, 0) . (30) 



Clearly, these equations have the structure similar to (7), 
but the inequalities (8) are satisfied only for sufficiently 
large A, that is, in the strong-coupling regime. Therefore, 
the continuum limit at A — > cannot be reached by the 
method described in Section 3. 

Let us, however, rewrite the equation (30) as 



i 



J2H^±^) + XS(x)Y(31) 



and introduce the "hopping parameter" 

1 



A + E£(±e M )' 



(32) 



Now the equation (30) in the form (31) looks like the 
equation for the free massive scalar propagator on the 



lattice with the mass m 2 = k^ 1 — 2D in lattice units. 
Note that in the weak-coupling limit A — > £ (±e M ) — ^ 1, 
m 2 — > and we approach the continuum limit. 

Let us now solve the equation (31) stochastically, as- 
suming that £ (x) is proportional to the stationary proba- 
bility distribution w (x) of some random process: £ (x) = 
cw(x), J2 W ( X ) = 1 ■ From (31) we get c = , A „ K „._ . The 

X 

equation (31) now looks as 



1-2Z5* 



w{x) = n^w{x±e tl ) + {l-2Dn) 8{x) . (33) 



Combining this equation with the definition (32), it is 
easy to show that k obeys the following self-consistency 
condition: 



(34) 



2D + Xw (0) ' 



The equation (33) has the form (7) without the nonlinear 
term and thus can be interpreted as the equation for the 
stationary probability distribution of the position of an 
ordinary bosonic random walk, defined by the following 
possible actions at each discrete time step: 

Move: With probability 2Dk move along the random 
unit lattice vector ±e„. 

Restart: With probability (1 — 2Dk) start again at the 
origin x = 0. 

This ensures that w (0) > and hence K never exceeds 
its critical value k c = (2D) . Therefore £ (a;) and w (x) 
can be expanded in powers of n. 

Thus we have defined a new expansion parameter k, 
which should obey the self-consistency equation (34), and 
obtained a well-defined convergent expansion, namely, 
the sum over all paths on the lattice with the weight k l , 
where L is the length of the path. Note that the quantity 
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£ (±e M ) in fact plays the role similar to the gluon conden- 
sate (which is expressed in terms of the mean plaquette in 
lattice theory) in non-Abelian gauge theory: one can ab- 
sorb all the divergences into the self-consistent definition 
of condensates [15, 16]. 

The final step in the construction of the nonlinear ran- 
dom process which solves the equations (30) is the solu- 
tion of the self-consistency equation (34). One possible 
solution is to use the iterations 



1 



2D + Xw (0; Ki) ' 



(35) 



Here w (0; k) is the return probability of a bosonic ran- 
dom walk with hopping parameter k. In practice, one 
should simulate the bosonic random walk at fixed k = Ki 
for some number T of discrete time steps, and then es- 
timate w (0; Ki) as w (0; «») w t (0) /T, where t (0) is the 
number of discrete time steps spent at x = 0. From (35) 
one then gets and the process is repeated until the 

value of k stabilizes with sufficient numerical precision. 
We call such algorithm "Algorithm A" . One can also con- 
sider an ultimate case, for which the return probability 
is updated and estimated as t(0) ft every time the point 
x = is reached. Now t is the time from the start of 
the random process and t (0) is the number of time steps 
spent at x = 0. Such algorithm will be called "Algorithm 
B". 

Mathematically, such random processes are not 
Markov processes, since the transition probabilities at 
each next step depend (via Ki) on the behavior of the 
process at all previous time steps. Stationary probability 
distributions of such processes obey nonlinear equations 
(such as (30)) [8], and and hence they are also called 
nonlinear random processes. 

As an interesting side remark, let us discuss such a 
theory at finite temperature, which is described by a 
bosonic random walk on the cylinder. Clearly, an ordi- 
nary bosonic random walk does not feel this compactifica- 
tion of space, and its stationary probability distribution 
is just a periodic linear combination of the correspond- 
ing distribution in infinite space. Such behavior cannot 
lead to any nonlinear finite-temperature effects such as 
phase transitions. On the other hand, if the parameters 
of the random walk depend on the return probability, as 
in (34), there is a nonlinear feedback mechanism since in 
the compactificd space the returns are more likely. Thus 
finite temperature indeed affects the local behavior of the 
random walker with memory and might lead to interest- 
ing critical phenomena. 

In order to illustrate such a stochastic solution of 
the equations (30), we consider the case D = 2. In 
two dimensions the model (29) is asymptotically free, 
and one can introduce the lattice spacing by fixing the 
value of mass in physical units (we set m p hy s = 1) : 
m p hys a (A) = m /att (A) = \J 'k -1 (A) - 2D. The process 
of convergence of the lattice spacing to its exact value 
is illustrated on Fig. 11 for both the algorithms "A" 
and "B". For algorithm "A" we have used T = 5 ■ 10 5 . 
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FIG. 11: The process of convergence of the random process 
with memory which solves the Schwinger-Dyson equations 
(30) for D — 2. The quantity plotted is the estimate of the 
lattice spacing a (A, t) after t discrete time steps, with the 
exact result o (A, t — I oo) subtracted. 
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FIG. 12: Numerical estimates of the lattice spacing as a func- 
tion of the coupling constant for the two-dimensional large- iV 
O (N) sigma model, compared with the exact result. The es- 
timates were obtained using both algorithms "A" and "B" 
with different number of time steps. 



The algorithm "A" converges much faster than the algo- 
rithm "B" . The values of lattice spacing obtained using 
both algorithms are compared with the exact solution on 
Fig. 12. In agreement with asymptotic freedom, lattice 
spacing quickly decreases with A. Again, algorithm "A" 
yields more precise results in the same number of time 
steps. 



Discussion and conclusions 

In this paper we have presented numerical strategies 
for the stochastic summation and re-summation of pcr- 
turbativc expansions in large- N quantum field theories. 
Our basic approach was to interpret the Schwinger-Dyson 
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equations as the equations for the stationary probability 
distribution of some random process. Since Schwinger- 
Dyson equations in such theories are nonlinear equa- 
tions, we had to use so-called nonlinear random pro- 
cesses, rather than ordinary Markov processes whose sta- 
tionary probability distributions always obey linear equa- 
tions. 

It is interesting to note that since the configuration 
spaces of random processes described in this paper are 
discrete, their numerical implementation require floating- 
point operations only for the random choice of actions. 
Thus such algorithms can be potentially much faster than 
the standard Monte-Carlo simulations based on floating- 
point arithmetic, and can be advantageous for machines 
based on CPUs. 

Our final goal is to extend the presented approach to 
non-Abelian lattice gauge theories. However, in this case 
direct stochastic interpretation of Schwinger-Dyson equa- 
tions is only possible at strong coupling, while the con- 
tinuum limit of such theories corresponds to the weak- 
coupling limit. In Section 4, we have discussed a way 
to access the weak-coupling limit, which, however, was 
implemented numerically only for O (N) sigma-model at 
large N. The basic idea is to absorb the divergences 
into a self-consistent redefinition of the expansion pa- 
rameter and solve the self-consistency conditions using 
random processes with memory. In some sense, O (N) 
sigma-model can be thought of as the bosonic random 
walk in its own condensate, and the approach to the 
self-consistent value of mass gap (see Fig. 11) - as a 
renormalization-group flow. 

For non-Abelian gauge theories the redefined expan- 
sion parameters can emerge as the lagrange multipli- 
ers for the "zigzag symmetry" of the QCD string and 
should also satisfy some self-consistency conditions [16]. 
Zigzag symmetry means that when one adds a line which 
is passed forward and backward to the boundary of the 
fluctuating string, the amplitudes should not change. In 
lattice gauge theory, this condition is equivalent to the 
unitarity of the link variables U {x, fi), which is similar to 
the condition \n (x) \ = 1 in O (N) sigma-model. These 
redefined parameters can be also related to the gluon 
condensate [16]. By analogy with the sigma-model, one 
can think that non-Abelian gauge theories are similar 
to strings moving in some self-consistent condensates. 
Such a picture is also close to the idea of holographic 
AdS/CFT duality for non-Abelian gauge theories, where 
the dual string lives in some self-consistent gravitational 
background, and the parameters of this background can 
be related to gluon condensates in gauge theory [11]. In 
fact, the requirement that the metric of the holographic 
background approaches that of the AdS space-time en- 
sures the zigzag symmetry of the strings which end on 
the AdS boundary [11]. 

In view of these qualitative considerations, our hope 
is that the loop equations in non-Abelian gauge theories 
can be solved stochastically by a random process similar 
to the one which was devised for the Weingarten model 



of random surfaces (see Subsection 3.3), but with some 
self-consistent choice of parameters, which might be im- 
plemented as the "memory" in the random process. 

Among other possible applications of the presented 
method one can think of the solution of Schwinger- 
Dyson equations in continuum gauge theories, combined 
with the Renormalization Group methods [10], numerical 
analysis of quantum gravity models described by various 
matrix models, and numerical solution of hydrodynami- 
cal equations [23]. 

It should be noted here that several attempts at the 
stochastic solution of the loop equations in large- N gauge 
theories have been already described in the literature 
quite a long time ago [22]. These algorithms were, in 
essence, based on the so-called branching random pro- 
cesses, so that the Wilson loop W (C) is proportional to 
the probability of transition from the initial loop config- 
uration C to the empty configuration with no loops. In 
particular, in contrast to the algorithm described in Sub- 
section 3.3, where one of the basic steps is to join loops, 
in the algorithms described in [22] the basic step was to 
split a self-intersecting loop into two loops. As a result, 
these algorithms did not implement the importance sam- 
pling and were not able to produce any sensible results for 
the four-dimensional gauge theory. Generally, branching 
random processes similar to those considered in [22] can 
be obtained from the "recursive" nonlinear random pro- 
cess described in this paper by time reversal. However, 
since such processes do not satisfy any detailed balance 
condition, they are not invariant under this operation, 
and lead to very different numerical algorithms. 
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Appendix A: Stochastic solution of nonlinear 
equations with coefficients of arbitrary sign 

The random process described in Section 2 was de- 
vised under the assumption that the coefficients p c (x) , 
Pe(x\y\) and pj (x\y±, 1/2) are real and positive for any 
x, yi and 2/2 • l n this Appendix we show how the solu- 
tion of the equation (7) with arbitrary signs or complex 
phases on the r.h.s. can be reduced to the solution of 
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another equation of the form (7) with all positive coef- 
ficients, provided the inequalities (8) are satisfied. We 
begin by discussing the case of real but non-positive co- 
efficients in details, and finally sketch the extension to 
complex- valued coefficients. 

To this end, let us represent the sign-alternating coef- 
ficients in (7) as: 

Pc (x) = P c +) (x) - p[r ] (x) 
Pe (x\y) = p { e +) (x\y) - p{- } (x\y) 
Pj (x\y) = Pj (x\y 1 ,y 2 ) - p\ r) (x\y 1 ,y 2 ) , (Al) 



ft is now easy to check that the difference 

w (x) = (x) - (x) (A4) 

satisfies the equation (7). On the other hand, the equa- 
tions (A3) again have the form of (7) with all positive 
coefficients, but with the configuration space X' being 
the direct product X ® 1 2 . In other words, each variable 
x now in addition carries the "sign", + or — , which can 
be written as {x, +} or {x, — }. Basing on the results pre- 
sented in Section 2, one can devise the random process 
which solves these equations: 

Create: With probability pc (x) create new clement 
{x, ±} £ X' and push it to the stack. 

Evolve: With probability p e + ^ (x\y) pop the clement 
{y, ±} from the stack and push the element {x, ±} 
to the stack, with probability pi ^ (x\y) do the same 
but flip the sign of y. 

Join: With probability p^ {x\yx,yz) consecutively pop 
two elements {j/i,si}, {y 2 ,s 2 } from the stack and 
push a single element {x, sis 2 } to the stack. That 
is, two pluses or two minuses associated with y's 
give {x, +}, but one plus and one minus give {x, — }. 



where p c (x), p e (x\y) and p^ (x\yi,y2) are all posi- 
tive and also obey the following inequality: 

J2 Y,Pc } (*) +Pi s) (x\vi)+pP {x\yx,y 2 ) < 1 (A2) 

for any yi , y 2 ■ Obviously, these inequalities can be satis- 
fied if the inequalities (8) are satisfied. Let us now intro- 
duce two functions (x) and if/ - ' (x), which satisfy 
the following equations: 



(A3) 



With probability p ( ~ > (x\yi,y 2 ) do the same, but 
flip the resulting sign, that is, push the element 
{x, — sis 2 } to the stack. 

Restart: Otherwise, empty the stack and push a single 
element {x, ±} € X' into it, with probability pro- 
portional to pc (x). 

As in Section 2, (x) and u> ( - _ - ) (x) are proportional to 
probabilities to find the elements {x, +} or {x, — } on the 
top of the stack, provided there is more than one element 
in it. 

The extension of this construction to complex- valued 
coefficients is quite straightforward. We represent the co- 

efficients in (7) as p c (x) = J d9 p c (x,8) exp {id), where 



p c (x, 9) is real and positive, and similarly for the other 
coefficients. The configuration space X' becomes the di- 
rect product X (g) S 1 , where S 1 is the unit circle in the 
complex plane. The stack now contains the pairs {x,8}, 
with 8 G [0, 2ir] being the complex phase. The func- 

tion w (x) is estimated as w (x) = J ddw (x, 8) exp {iff), 



where w (x, 8) is the probability to find the element {x, 9} 
at the top of the stack. In the random process, new el- 



™ (+) (x) = P w (x) + J2pi +) ( x \y) wi+} 0) + (%) wi ~ ] (y) + 

y y 

^2 p< j +) {x\yi,y 2 )w {+) (yi)w {+) (y 2 ) + ^2 {x\yi,y 2 )w { ~ ] (yi)u> (_) (y 2 ) + 
n,V2 yi,V2 

+ J! pf {x\yi,y 2 ) w {+) (j/i) w { ~ ] (2/2) + ^2 p i 2/2) (y x ) w i+) (y 2 ) 

yi,y2 yi,V2 

(X) = p(~) (X) + J2Pe +) ( X \v) Wi ~ } (V) + Y. P e~ ] ^ wi + ) (») + 

y y 
^2 Wfi'Sfc) w {+) (yi)w ( ~ ) (y 2 ) + ^2 p\ +) {x\yi,y2)w ( -~ ) (yi)w (+) (y 2 ) + 

1/1,1/2 Vl,V2 

+ ^2 p i ( x \y^y*) wi+) (2/1) (2/2) + ^2 V Y ] ^) w { ~ ] (yi) (y 2 ) 
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ements are created with probability distribution propor- 
tional to p c (x,6), and in the "Evolve" and the "Join" 
actions the phases of the elements and the coefficients in 
(7) are added modulo 2ir, similarly to signs. 

Note that this solution does not in general have the 
property of importance sampling. Indeed, one can have 
some x for which w (x) is numerically very small, but the 



random process can spend an almost equal large amount 
of time in the states {x, +} and {x, — }, so that numeri- 
cally large (x) and (x) nearly cancel. Whether 
this occurs or not depends on the particular system of 
equations and on the particular unknown variables, but 
potentially this feature can make numerical simulations 
less efficient. 
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